Use Pearson residuals for glmmTMB/MixMod heteroscedasticity plot - #460
Use Pearson residuals for glmmTMB/MixMod heteroscedasticity plot#460jbogomolovas2 wants to merge 1 commit into
Conversation
There was a problem hiding this comment.
Pull request overview
Updates see’s heteroscedasticity plotting for glmmTMB/MixMod models to use Pearson residuals (matching the intended variance-standardization behavior discussed upstream in easystats/performance), and documents/releases the fix.
Changes:
- Prefer
residuals(type = "pearson")forglmmTMB/MixModinplot.see_check_heteroscedasticity(), with a fallback to the prior scalar-scaling approach. - Add a NEWS entry describing the bug fix.
- Bump package version and add a contributor entry in
DESCRIPTION.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 2 comments.
| File | Description |
|---|---|
| R/plot.check_heteroscedasticity.R | Switch heteroscedasticity plot residuals for glmmTMB/MixMod to Pearson residuals with fallback logic. |
| NEWS.md | Document the bug fix in the changelog. |
| DESCRIPTION | Increment version and add contributor metadata. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| sig <- if (faminfo$is_mixed) { | ||
| sqrt(insight::get_variance_residual(model)) | ||
| } else { | ||
| .sigma_glmmTMB_nonmixed(model, faminfo) | ||
| } | ||
| stats::residuals(model, type = "response") / sig |
There was a problem hiding this comment.
@jbogomolovas2 I think this is the same as for performance, which should be fixed. Else, everything looks good to me.
There was a problem hiding this comment.
SO, basically, just change it to "betadisp".
| r_pearson <- tryCatch( | ||
| stats::residuals(model, type = "pearson"), | ||
| error = function(e) NULL | ||
| ) | ||
| if (is.null(r_pearson) || all(is.na(r_pearson))) { |
There was a problem hiding this comment.
I'll add visual tests later
| sig <- if (faminfo$is_mixed) { | ||
| sqrt(insight::get_variance_residual(model)) | ||
| } else { | ||
| .sigma_glmmTMB_nonmixed(model, faminfo) | ||
| } | ||
| stats::residuals(model, type = "response") / sig |
There was a problem hiding this comment.
@jbogomolovas2 I think this is the same as for performance, which should be fixed. Else, everything looks good to me.
Companion to easystats/performance#927, which @strengejacke asked me to open.
The problem
plot.see_check_heteroscedasticity()carries its own copy of the branch fixedin performance#927, including a private copy of
.sigma_glmmTMB_nonmixed(). Soplot(check_heteroscedasticity(m))is affected the same way: forglmmTMBandMixModmodels the residuals are divided by a single scalar, which is onlycorrect when
V()does not depend onmu. For non-mixed binomial and Poissonmodels that scalar is
1, so no standardization happens and the panel bows forcorrectly specified models.
The full write-up, including the per-family table showing this affects nine of
eleven glmmTMB families rather than just binomial, is in performance#927.
The change
Same fix, ported. Try
residuals(type = "pearson"), keep the existingexpression as a fallback.
seehas no.safe()helper, so this usestryCatch()directly, matching the surrounding code.Effect
Correctly specified glmmTMB binomial fit (n = 600, size = 20). Left: patched.
Right: current behaviour, where the trend line humps in the middle purely
because binomial spread is largest near p = 0.5.
Also worth flagging,
.sigma_glmmTMB_nonmixed()in R/plot.check_heteroscedasticity.R readsglmmTMB renamed that parameter to
betadisp(glmmTMB commit 472ccbe0);performance's copy of this helper was updated, this one was not.
[does not partial-match, so the lookup returns NA,exp(NA)is NA, anddividing by it makes the whole panel NA. Nothing errors, so the surrounding
tryCatch does not catch it and the plot comes out silently blank.
Affects gaussian and Gamma glmmTMB fits. binomial/poisson/truncated_poisson
short-circuit to 1 before reaching that line.
Confirmed on glmmTMB 1.1.15:
Fix is the rename, in all three places in the switch(). Happy to PR it, but
noting it separately from # since it is a different failure mode.
Notes
performance#927.